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MODIFICATIONS  AND  IMPROVEMENTS  IN  A 
STRUCTURAL  OPTIMIZATION  SCHEME  BASED 
ON  AN  OPTIMALITY  CRITERION 


1 . INTRODUCTION 

While  recent  years  have  seen  continuing  increases  in  the  variety 
and  scope  of  applications  of  structural  optimization  technology  within 
the  aerospace  industry,  there  are  still  some  significant  problems  to 
be  overcome  before  this  technology  can  be  routinely  applied  wherever 
it  is  needed  in  all  stages  of  aerospace  vehicle  design.  Among  these 
problems  are  the  need  for  efficient  treatment  of  large  numbers  of 
design  variables  (in  the  thousands,  for  example)  with  many  different 
constraints.  Indeed,  it  may  well  be  that  the  number  of  different 
design  requirements  with  which  industry  is  faced  is  growing  more 
rapidly  than  the  number  of  constraints  that  is  being  incorporated  in 
the  most  recent  structural  optimization  schemes. 

In  an  effort  to  aid  in  the  resolution  of  these  problems,  Nielsen 
Engineering  & Research,  Inc.  (NEAR)  has  underta)cen  a research  program 
under  the  sponsorship  (Contract  F49620-77-C-0055)  of  the  Air  Force 
Office  of  Scientific  Research  (AFOSR) . Of  specific  interest  is  the 
extension  of  an  optimality-criterion  algorithm  (refs.  1-3)  to  large 
problems  involving  multiple  constraints,  where  the  constraints  include 
both  strength  and  stiffness  requirements.  Before  proceeding  with  this 
algorithm,  however,  it  was  considered  advisable  to  conduct  some  effi- 
ciency comparisons  in  order  to  identify  other  algorithms  with  equal 
or  perhaps  superior  qualities.  Other  tasks  that  were  to  be  undertaken 
during  the  first  year  involved  a study  of  the  accuracy  of  computing 
flutter-speed  or  flutter-eigenvalue  gradients  with  fixed  modes  defining 
generalized  coordinates  and  possibly  the  incorporation  of  multiple 
equality  behavioral  constraints.  The  sections  that  follow  describe 
the  first  year's  accomplishments  in  more  detail,  the  conclusions  to  be 
drawn  from  this  work,  and  the  work  planned  for  the  next  year.  One 
archive  publication,  covering  a portion  of  the  year's  accomplishments, 
will  appear  in  1978.  A preprint  of  this  paper  is  appended. 
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2.  ALGORITHMS  BASED  ON  OPTIMALITY  CRITERIA 

2.1  Optimality  Criteria  Versus  Mathematical  Programming 

Optimization  algorithms  are  more  or  less  arbitrarily  classed 
according  to  the  philosophy  underlying  their  derivation.  Mathematical- 
programming  (MP)  algorithms  employ  a relatively  sophisticated  set  of 
calculations  at  each  design  point  in  order  to  produce  a design  change 
that  will  simultaneously  improve  the  merit  function  (i.e.,  reduce  the 
total  weight)  without  violating  any  of  the  constraints.  Optimality- 
criterion  (OC)  algorithms  are  derived,  generally  heuristically , from 
conditions  that  must  be  satisfied  at  the  optimum  design.  MP  algorithms 
can  generally  be  proven  to  converge,  while  OC  algorithms  cannot.  On 
the  other  hand,  OC  algorithms  are  less  involved  than  MP  algorithms  from 
a computational  standpoint.  However,  the  fundamental  reason  for 
interest  in  OC  algorithms  is  to  be  found  in  their  potential  for  applica- 
tion to  very  large  problems,  based  originally  on  the  success  of  the 
stress-ratio  or  f ully-stressed-design  algorithm  in  treating  stress 
constraints  on  members  numbering  in  the  thousands. 


2.2  Optimality  Criterion  and  Recursion  Relation  for  an  Equality 
Behavioral  Constraint 

OC  algorithms  come  in  many  forms.  To  fix  ideas,  let  us  consider 
a single  behavioral  constraint  on  a structure  to  be  optimized.  The 
structure  is  characterized  by  N design  variables  t^,  and  its  mass  m 
to  be  minimized  is  assumed  to  be  a linear  function  of  these  variables; 


m(t^) 


m + m(t . ) 
o 1 


m + 
o 


N 


i=l 


m.  t . 
1 1 


(1) 


The  assumption  of  linearity  is  not  necessary,  but  it  covers  a broad 

class  of  finite-element  models  where  the  design  variables  are  plate 

thicknesses,  bar  areas,  or  the  like.  In  addition,  there  is  some  mass 

m that  is  not  available  for  optimization,  which  may  represent 
o 

fasteners,  lightly  loaded  structure,  fuel,  etc.  The  behavioral 
constraint  is  written  simply  as 

c(t^)  = 0 


(2) 


and  there  may  also  be  minimum-gage  constraints  of  the  form 


1 1 min 


(3) 


One  simple  form  of  an  optimality  criterion  can  be  obtained  from  varia- 
tions of  the  merit  function 


J = m + Xc 


(4: 


In  particular,  the  requirement  of  vanishing  variations  of  J with 
respect  to  the  active  design  variables  leads  to 


1 9c 


m.  9t . 
1 1 


1 

X 


V i e A 


(5) 


Here  A denotes  the  set  of  design  variables  that  are  active — i.e.,  not 
at  their  minimum  values--at  the  optimum.  A number  of  recursion  rela- 
tions have  been  introduced  for  iteratively  resizing  the  structure  in 
order  to  arrive  at  a design  that  satisfies  equation  (5) . (A  very 
revealing  discussion  of  the  relationships  among  many  of  these  rela- 
tions may  be  found  in  ref.  4.)  One  that  has  been  widely  used  is  due 
to  Kiusalaas  (ref.  5) : 


where 


t'’ 

1 11 


(1 


a.  ) X 


(6) 

(7) 


Here  v is  the  iteration  number,  and  a is  a parameter  that  ranges 
value  from  -1  to  +1.  it  may  or  may  not  vary  with  each  iteration. 


Note  that  as  the  optimum  is  approached,  X 
no  matter  what  value  a'^  takes. 


9c 


i 


-♦  -1 , and  C . 

1 


in 
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2.3  Two  Resizing  Algorithms 

There  are  in  turn  a number  of  ways  of  devising  resizing  algorithms 
based  on  equations  6 and  7.  One  procedure,  from  references  1-3,  deter- 
mines a'^  so  that  a given  reduction  in  mass  is  achieved,  while  the 


r 


multiplier  A is  determined  from  the  requirement  that  the  constraint 

v 

value  c be  always  driven  to  zero.  To  first  order,  the  change  in  the 
constraint  value  is  given  by 


^v+ 1 V 

c - c 


_ r [ f v+1  V 

- l3vr)  H 


With  Am 


N 

V V / 1 

~ L ”’•1  • equations  (5) -(8)  can  be  manipulated  to 

i=l  ^ ^ 1 


yield,  for  c^^^  = 0, 


V V V . V 

Am  + c 8^/62 

77^7277^ 


8^  + cV(a^  - 1) 

A ^ 

ft'' 


where 


ft''  - T f 3c  . I"*  H ?!  If  3er1  ,v 

^1  - J,  at-  ^ "2  = in-  FtT  ^i- 

1=1  ^ 1 ■'  1=1  1 1 


The  mass  change 


Am  is  specified  as 


Am  = -Km 


where  K is  chosen  from  a table  of  user-selected  values  arranged  in 
descending  order.  The  largest  value  is  chosen  first,  and  a^,  X^,  and 
the  cV  are  calculated  for  the  active  set  A.  If  any  active  design 
variable  is  changed  by  more  than  25%,  the  next  lower  value  of  K is 
chosen,  and  the  process  is  repeated  until  no  design-variable  change 
exceeds  the  25%  limit.  Next,  the  new  values  of  the  active  design 
variables  are  compared  to  their  minimum  gages.  If 

the  next  lower  value  for  K is  selected,  and  the  redesign  process  is 
repeated.  If  tV^^  < (t.)  for  the  smallest  value  of  K,  then  for 

that  design  variable  tV^^  = (t.)  . , and  it  is  relegated  to  the 

11  min 


passive  set.  A final  pass  through  the  redesign  stop  is  now  required. 


r 
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since  equations  (9)  and  (iO)  must  be  altered 

accordingly  to  reflect  revised  expressions  for  o'"  ^ - c^  and  Am'^  (see 
reference  1 for  details).  If  the  active-passive  identities  are 
unchanged,  then  a new  iteration  is  begun.  Convergence  is  checked 
in  two  ways.  First,  if  K is  at  its  minimum  value  and  an  active 
design  variable  is  calculated  to  change  by  more  than  25%,  the  design 
is  declared  final,  subject  to  a check  on  the  active-passive  identities. 
This  test  is  essentially  a "diminishing  returns"  criterion,  since  a 
mass  reduction  of  a given  amount  requires  larger  and  larger  chanyt-s 
in  the  design  variables  as  the  optimum  is  approached.  The  second 
test  involves  the  proximity  of  the  redesign  factors  cV  to  unity: 

\cl  - 11  < e (12) 


Here  e is  a user-supplied  convergence  parameter.  If  this  test  is 
satisfied,  the  design  is  declared  final,  again  subject  to  a check  on 
the  active-passive  identities.  This  check  is  to  ensure  that  there  are 
no  passive  design  variables  that  should  be  reintroduced  into  the  active 
set.  It  is  based  on  the  Kuhn-Tucker  optimality  conditions  and  is 
described  in  detail  in  references  1 and  . A flow  diagram  for  this 
algorithm  is  given  in  Appendix  A. 


A variant  on  the  above  algorithm  involves  simply  scheduling  a, 

rather  than  computing  it  for  a scheduled  set  of  mass  reductions 

(refs.  6 and  7).  Also,  the  minimum-gage  constraints  are  handled 

differently.  At  each  iteration  following  the  one  where  a design 

variable  reaches  its  minimum  value,  the  factor  cY  for  this  variable 

1 

is  still  evaluated.  If  > 1,  the  variable  is  reintroduced  immediately 
into  the  active  set.  (Here  also  special  steps  must  be  taken  in  the 
calculation  of  to  account  for  design-variable  changes  associated 
with  entry  to  or  exit  from  the  passive  set.  Details  are  in  references 
6 and  7.)  In  this  algorithm,  then,  active-passive  identities  are 
continually  checked,  but  there  is  very  little  computation  associated 
with  determining  This  parameter  is  scheduled  by  the  following 

formula,  which  replaces  equation  (9) : 

V v-1  (0)  , . v-1 

a = a Cl  =01  ( a ) 

X X 


(13) 


Typically,  will  be  chosen  to  be  near  unity  - say,  0.9.  Then  at 

each  iteration  the  current  value  of  a is  obtained  by  multiplying  the 
previous  value  by  a factor  which  is  slightly  less  than  unity  - 
say,  0.95.  The  reduction  of  a at  each  step  is  equivalent  to  increas- 
ing step  sizes  as  the  optimum  is  approached,  in  order  to  accelerate 
convergence.  Convergence  is  evaluated  by  checking  the  uniformity  of 
the  weighted  derivatives,  according  to  equation  (5) : 


1 + X 


3c 


m. 

1 


at. 

1 


< e< 


¥ i e A 


(14) 


The  convergence  parameter  e is  user-supplied.  A flow  diagram  for 
this  algorithm  is  also  given  in  Appendix  A. 


2.4  An  Alternate  Approach  to  an  Optimality  Criterion 

The  form  of  the  optimality  criterion  derived  in  subsection  2.2 
is  very  general,  and  this  generality  is  easily  extended  to  multiple 
constraints  (ref.  5) - However,  it  will  prove  instructive  to  consider 
at  least  one  other  approach,  which  involves  using  the  equations 
governing  the  structure  for  the  particular  constraint  being  considered. 
In  an  effort  to  retain  some  generality,  let  the  equations  be  written 
as 

lB(n,t^)]{q}  = {0}  (15) 


Here  {q}  is  a column  of  unknowns  which  may  be  discrete  displacements 
or  modal  coordinates,  and  the  coefficients  of  the  equations  are  given 
in  [B] , which  is  a function  of  the  design  variables  t^  and  possibly 
an  eigenvalue  as  well.  In  the  case  of  an  aeroelastic  constraint, 
such  as  a fixed  flutter  speed,  then  equation  (15)  will  be  complex; 
for  a constraint  on  a natural  frequency,  the  equations  will  be  real. 

The  merit  function  J can  now  be  written  as 

J = m + Re(LpJ  [Bl{q})  (16) 

Here  the  scalar  multiplier  X is  replaced  by  a set  of  multipliers  {p}, 
and  the  real  part  of  the  triple  matrix  product  must  be  taken  if 
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equations  (15)  are  complex  (ref.  8).  At  the  optimum,  the  merit 
function  J is  stationary  with  respect  to  independent  variations  of 
the  t^,  {p},  and  {q}.  If  the  eigenvalue  f)  is  not  fixed  by  the 
constraint,  then  variations  of  it  must  be  considered  also.  Applying 
these  conditions  yields 


[B]{q}  = [Oj 

(17) 

[B]  = (_0j  , or  [B]'^{p}  = {0} 

(18) 

(19) 

mi  + Re(|_pj[3^3^q})  = 0 

(20) 

Equation  (20)  is  the  optimality  criterion.  The  term  involving  the 
triple  matrix  product  resembles  an  energy  density,  so  this  could  be 
referred  to  as  an  "energy-density"  form  of  an  optimality  criterion. 
(Specific  expressions  for  various  types  of  constraints  will  be  pre- 
sented in  Section  3 below.)  Writing  Re  ([p|  (q)) 

permits  equation  (20)  to  be  recast  as  ^ 

(ev)i  = -1/  V i e A (21) 


As  noted  here,  in  the  presence  of  side  (e.g.,  minimum-gage)  constraints 
only  the  indexes  for  the  active  design  variables  are  to  be  considered. 


2.5  "Energy-density"  Recursion  Relation  and  Redesign  Algorithm 


The  number  of  recursion  relations  that  can  be  devised  to  satisfy 
equation  (21)  is  limited  only  by  the  designer’s  imagination.  One 
choice  is  a variation  of  the  recursion  relations  developed  in  refer- 
ences 9 and  10: 


t 


V 

i 


where 


= 

1 


(e  )Y 

V 1 


<«av) 


(1  + 


(22) 


and 
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(e 


)^  = — y 

av  NA  . ^ 


NA 

leA 


(23) 


This  relation  differs  from  those  in  references  9 and  10  in  several 
important  respects.  In  reference  9,  which  is  concerned  with  altering 
a strength  design  to  meet  a flutter-speed  constraint,  only  the  strain- 
energy  density  is  used  to  define  the  (e^) . , rather  than  the  complete, 
or  exact,  expression  represented  by  equation  (20) . Also,  the  denom- 
inator in  the  "energy-density"  ratio  is  obtained  by  considering  only 
those  elements  whose  (e^) ^ exceed  the  average  (e^^) . This  has  the 
effect  of  permitting  only  increases  in  design  variables,  so  that  a 
particular  design  variable  that  is  too  large  for  strength  require- 
ments in  the  redesigned  structure  cannot  be  reduced.  This  limitation 
is  removed  in  reference  10,  but  "energy-densities"  are  defined  as 


'Vi  * 


LpJ 


9B 

ati 


{q} 


which  is  also  different  from  that  required 


to  satisfy  equation  (20)  . Hence  neither  of  these  two  procedures  is 
capable  of  converging  to  the  "exact"  optimum,  whereas  a redesign 
algorithm  based  on  equations  (6),  (22),  and  (23)  will  do  so. 


In  equation  (23) , the  average  (e  ) 

3 V 


is  determined  by  averaging 
the  energy  densities  only  for  those  design  variables  in  the  active 
set.  (NA  is  the  number  of  active  design  variables.)  In  equation  (22), 
e^^  is  typically  less  than  unity  (for  example,  0.5)  and  02  is  greater 


than  unity  (say,  2.0). 


Since  the  individual  (e^)j^  n>ay  differ  in  sign 


from  the  absolute  value  of  the  "energy-density"  ratio  is  required 

to  avoid  numerical  problems.  This  ratio  also  implies  that  a design 
variable  should  be  decreased  only  if  its  "energy  density"  is  less 
than  the  average.  In  the  factor  involving  the  current  constraint 
value  o'*,  c^  > 0 represents  an  infeasible  value,  requiring  an  increase 
in  the  design  variables. 

A redesign  algorithm  based  on  equations  (6) , (22) , and  (23)  is 
quite  simple. 


Minimum-gage  constraints  are  treated  by  simply  relega- 
ting a design  variable  to  the  passive  set  whenever  tV^^  i ^^i^min* 
each  step,  all  cV  are  calculated,  so  whenever  a particular  cV  is 
greater  than  unity  for  a passive  design  variable,  that  variable  is 
reintroduced  into  the  active  set  and  the  redesign  procedure  is  invoked 


-9- 


again.  This  approach  to  the  identification  of  the  active  and  passive 
sets  is  identical  to  that  followed  in  the  algorithm  of  references  6 
and  7.  Convergence  is  evaluated  by  equation  (12).  A flow  diagram 
for  this  "energy-density"  algorithm  is  given  in  Appendix  A. 

3.  CONSTRAINT  EVALUATIONS  AND  DERIVATIVE  OR  "ENERGY-DENSITY" 
CALCULATIONS 

The  algorithms  discussed  in  Section  2 above  were  coded  as  separate 
routines  that  require  constraint  evaluations  and  calculations  of  the 
derivatives  or  the  "energy-densities".  In  this  section,  specific 
forms  for  the  constraints  considered — displacements,  natural  frequen- 
cies, and  flutter  speeds — are  presented.  Also,  the  equivalence  of  the 
two  forms  of  optimality  criteria — equation  (21)  or  equation  (5) --is 
discussed. 

3.1  Displacement  Constraint 

The  displacements  of  a structure  loaded  by  a set  of  forces  {F} 
are  found  by  solving 

[K]{u}  = {F}  (24) 

where  (K]  is  the  discrete  stiffness  matrix,  {f}  is  a column  of  nodal 

applied  forces,  and  {u}  is  a column  of  nodal  displacements.  One  of 

V ^r 

these  displacements,  u , is  to  be  constrained,  so  that  c = j- — t 1. 

^ ^“r'des 

In  the  derivative  calculation,  the  dummy-load  method  is  used.  This 

involves  calculating  the  displacements  {u  } due  to  a dummy  load  set 

{F^^^ },  where  F^^^  = 6^^: 

[kKu^"^^  } = {F^^^  } (25) 

It  is  now  assumed  that  the  stiffness  matrix  [K]  is  a linear  function 
of  the  design  variables  and  can  be  written  as 

N 

IK]  = [K  ] + I t.  IK.] 

° i=l  ^ ^ 


(26) 
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The  justification  for  this  assumption  follows  the  same  line  of 
reasoning  used  in  justifying  equation  (1).  The  constraint  derivative 
can  then  be  shown  to  be  (see,  for  example,  ref.  5) 


3c 

■5t:  ~ 

1 


(U 


r des 


r des 


(27) 


The  first  form  of  the  optimality  criterion,  equation  (5)  , becomes 


m.  (u  ) , 

1 r'des 


(_uj  [K^]  {u^^^  } 


1 

X' 
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(28) 


While  the  form  of  the  constraint  equation  in  this  case  does  not 
correspond  to  that  used  in  developing  the  second  form  of  the  optimality 
criterion,  equation  (21)  , it  is  nevertheless  possible  to  ma)ce  use  of 
the  recursion  relations,  equations  (6) , (22) , and  (23) , by  replacing 

(e  ) . with  . This  is  reminiscent  of  the  recursion  relation 

V 1 o t . 

1 

finally  chosen  for  FASTOP  (ref.  11),  in  the  sense  that  it  can  be 
viewed  as  a heuristically  derived  recursion  relation  based  on 
satisfying  the  optimality  criterion  given  by  equation  (5) . 

Constraint  evaluation  in  this  case  is  straightforward.  The 
derivative  matrices  IK^^I  are  invariant,  as  a result  of  the  linearity 
assumption,  so  equation  (26)  is  used  to  update  the  stiffness  matrix 
at  each  step. 

The  displacement  u is  found  by  solving  equations  (24)  for  {u}. 

^ (r ) 

The  dummy- load  displacements  {u  } are  found  by  solving  equations 
(25)  , and  is  calculated  from  equation  (27)  . This  information  is 

then  passed  to  any  of  the  optimization  subroutines. 


3.2  Frequency  Constraint 

The  equation  of  motion  for  free  vibration  can  be  written  in 
discrete  form  as: 

( [K]  - 0)^  [M]  ) {u}  = {0}  (29) 

The  mass  matrix  is  also  assumed  to  be  a linear  function  of  the  design 
variables : 
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IM]  = [M^]  + I t.  [Ml  (30) 

w • « X X 


The  order  of  [M]  and  [K]  is  typically  large,  and  it  is  often  more 
efficient  to  rewrite  equation  (29)  in  terms  of  modal  coordinates, 
where  a set  of  natural  modes  of  the  initial  design  is  used  to  define 
generalized  coordinates  for  all  subsequent  designs.  (Although  this 
could  lead  to  inaccuracies  if  the  optimal  design  is  substantially 
different  from  the  initial  one,  periodic  updating  of  the  modes  or 
simply  using  a few  more  modes  of  the  initial  design  should  be 
sufficient  to  avoid  severe  problems.  This  is  an  open  question  that 
still  needs  to  be  answered.)  With  ti}>]  defined  as  the  modal  matrix, 
modal  coordinates  are  defined  as  {u}  = [4)]{q},  and  equation  (29) 
becomes 

([GK]  - [GM]) {q}  = {O},  (31) 


with 

[GK]  = [({.] 

IGM]  = [(}>] 

The  linearity  assumption  embodied  in  equations  (30)  and  (26)  also 
carries  over  to  the  generalized  mass  and  stiffness  matrices,  so  that 
derivative  generalized  matrices  [GM^l  and  [GK^]  may  be  defined  by 
obvious  analogy.  With  fixed  modes  defining  generalized  coordinates, 
these  are  also  invariant  and  of  much  smaller  order  than  their  discrete 
counterparts . 


If  0)  is  the  frequency  to  be  fixed,  and  (m)  , is  the  desired 
r , . x ces  I X ix 


(to) 


(0 


value  of  this  frequency,  then  c = 

3(0  " ~ ~r 

To  calculate  , equation  (31)  is  written  for  the  rth  mode,  differ- 


des  , 3c 

1,  and 

1 


(lo),  3(0 

des  r 

3t . • 
1 


to. 


3t. 
1 


entiated  with  respect  to  t.,  and  then  premultiplied  by  |_q ' J , the  rth 

^ 3(0 

eigenvector.  Solving  then  for  gives 


3t. 


3(0 

r 

9ti 


lGK.l{q'"M  - (o^[q^"^J  [GM.Hq'^M 


(r) 


(r)| 


(r) 


(r) 


T 


(33) 


-12- 


and  the  optimality  criterion  of  equation  (5)  reads 

iGK^llq^^^}  - [GM^]{q^^^}  = C,  V i £ A (34) 

where,  with  w = (to)  , , 

r des 

C = ^U'"’Jl=«l(.3‘"')  U5) 

The  left-hand  side  of  equation  (34)  can  be  termed  a "specific  Lagrangian 
density",  since  it  is  the  difference  between  the  pealc  strain  energy  and 
the  pea)c  Icinetic  energy  in  the  constrained  mode  per  unit  value  of  each 
active  design  variable.  Note  also  that  equation  (34)  is  homogeneous 
with  respect  to  the  eigenvector  {q  },  so  it  is  invariant  with  respect 

{ IT  ) 

to  the  normalization  of  {q  }. 

For  the  "energy-density"  form,  the  matrix  [B]  is  identified  with 
IGK)  - w [GM] , and  {p}  and  {q}  are  both  (q  since  [B]  is  symmetric 

(see  equations  (17)  and  (18)).  The  eigenvalue  is  fixed,  so  there 
is  no  "free"  eigenvalue  Q,  and  equation  (19)  does  not  apply.  Upon 
writing  out  equation  (21)  becomes 

Lq^^^J  [GK^liq^^^}  - V i e A (36) 

In  this  equation,  {q  ) can  be  normalized  to  give  an  arbitrary  value 
on  the  right-hand  side,  so  equations  (34)  and  (36)  are  in  fact 
equivalent. 

For  any  of  the  optimization  routines,  constraint  evaluation  is 

straightforward.  Equation  (31)  is  solved  for  w and  {q  },  and  then 

the  or  (e  ) . are  calculated.  The  mass  and  stiffness  matrices  are 

3ti  VI 

updated  with  the  current  values  of  the  design  variables  by  equations 
identical  in  form  to  equations  (30)  and  (26). 

3.3  Flutter-Speed  Constraint 

The  governing  equations  for  flutter,  written  in  modal  coordinates, 
have  the  form 


( [GM]  + [GA]  - fi[GK]){q}  = {0} 


(37) 
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Here  [GM]  and  IGK]  are  the  generalized  mass  and  stiffness  matrices 

defined  previously,  and  [GA]  is  the  generalized  aerodynamic  matrix. 

The  modal  coordinates  {q}  are  again  defined  by  a fixed  set  of  modes, 

so  [GA]  is  a complex  function  of  reduced  frequency  k and  Mach  number 

(The  altitude  is  assumed  given.)  The  eigenvalue  51  is  written 

2 

in  terms  of  the  frequency  w and  damping  parameter  g as  (1  + jg)/oj  , 
and  flutter  is  determined  by  that  combination  of  k and  which  gives 
g = 0 and  the  lowest  flutter  speed  U^,  which  is  obtained  from  k and  u>. 
A final  step  involves  varying  the  Mach  number  until  there  is  com- 
patibility among  the  altitude,  the  Mach  number,  and  the  flutter  speed. 

The  flutter  constraint  will  be  enforced  by  requiring  that  g = 0 
for  a given  combination  of  altitude,  Mach  number,  and  speed.  The 
derivative  is  calculated  in  the  manner  of  references  1 and  2: 


2-  = - gl^)  {2gR3  + 2X3  + ^ I^) 

OD 

- (2R3  - 2gl3  + ^ R4)  (ij  - + gR2)1/D 

00  -J 


where 


D = (2gR3  +213+^  1^)13  + (2R3  - 2gl3  + ^ R4)R3  (2 

^00  '^00 

rJ  = Re([pJ  IGM^I  {q}) 

R^  = Re(LpJ  [GK^J  {q}) 

R3  = Re([pJ  [GK]  {q})  {< 

R4  = 

and  I3,  1^,  I3,  and  are  the  corresponding  imaginary  parts.  It  is 
also  understood  that  {q}  is  the  eigenvector  for  the  critical  flutter 
mode,  and  {p}  is  the  eigenvector  for  the  adjoint  problem 


([GM]  + [GAl  - 5i[GKl){p}  = {0} 
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For  and  I^,  derivatives  of  the  elements  of  the  aerodynamic  matrix 
[GA]  are  needed  with  respect  to  reduced  frequency  k.  These  are 
obtained  approximately  by  calculating  [GA]  for  a band  of  reduced 
frequencies  in  the  range  of  interest  and  then  fitting  these  data 
with  polynomials  in  k.  The  fitting  procedure  used  here  is  the  one 
discussed  at  length  in  reference  12.  The  derivatives  are  then 
obtained  by  differentiating  the  polynomials.  For  a flutter  constraint, 
then,  the  o-'timality  criterion  is  obtained  from  equation  (5)  by 
identifying  c with  g: 


1 3q  _ 1 

m.  at.  X 


i 1 


f 
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with  given  by  equation  (38)  . 

For  the  ’’energy-density"  form,  [B]  is  identified  with  (GM)  + 
[GA]  - n[GK],  {p]  with  {p},  and  {q}  with  {q},  so  the  optimality 
criterion  in  this  form  is 


(45) 


^®v^i  " Re((_pJ  ( [GM^l  - niGK^l)  l^})  = -1,  V i e A 


(46) 


This  expression  is  much  simpler  than  equation  (45)  , so  it  appears  that 
an  algorithm  based  on  computing  (e„) . rather  than  would  be  more 

V 1 o t . 
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efficient,  since  less  computation  is  involved,  and  derivatives  of  the 
aerodynamic  matrix  need  not  be  calculated.  However,  this  is  only 
part  of  the  story,  since  the  convergence  characteristics  of  the 
algorithm  also  influence  the  efficiency. 


It  is  worth  noting  that  the  "energy-density"  formulation  provides 
an  additional  relation  when  a "free"  eigenvalue  is  involved.  In  this 
case,  the  "free"  eigenvalue  is  the  flutter  frequency,  which  is  uncon- 
strained. Thus  equation  (19)  , with  identified  with  u),  is  applicable. 
By  making  use  of  the  actual  expression  for  [B]  and  the  definitions  of 
equations  (40)- (43)  and  their  imaginary  counterparts,  equation  (19) 
can  be  manipulated  to  give  another  condition  to  be  satisfied  at  the 
optimum: 


2R3  - 2,13  t = 0 


(47) 


■i 

1 

\ 

i 
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If  this  relation  is  inserted  into  equations  (38)  and  (39),  equation 
(45)  is  greatly  simplified  and  reads 


i + ,1^)  . J , 


V i e A 


(48) 


Equation  (46)  , when  rewritten  with  the  definitions  of  equations  (40) 
and  (41)  and  their  imaginary  counterparts,  becomes 


1 t 2„i 

m.  1 

1 


^2 
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(49) 


Since  the  normalization  of  {p}  and  {q}  is  arbitrary,  it  is  always 
possible  to  renormalize  these  eigenvectors  so  that  the  right-hand 
side  of  equation  (49)  is  identical  to  the  right-hand  side  of  equation 
(48) . Hence  the  two  optimality  criteria  are  formally  equivalent. 
However,  it  must  be  emphasized  that  the  derivative  is  in  general 

d ti 

given  by  equation  (38)  ; the  simplified  expression  is  valid  only  at 
the  optimum. 

For  any  of  the  redesign  algorithms,  constraint  evaluation  is 
accomplished  by  solving  equation  (37).  Equation  (44)  is  solved  to 
obtain  {p}.  The  generalized  mass  and  stiffness  matrices  are  updated 
as  described  in  subsection  3.2.  Since  the  flutter  frequency  will 
vary  as  the  redesign  progresses,  the  reduced  frequency  k must  be 
updated  in  order  to  ensure  compatibility  between  its  value  and  the 
frequency  computed  from  fi.  This  is  most  easily  done  iteratively. 

The  value  of  k from  the  previous  design  is  used  initially  to  determine 
[GA] , and  a new  frequency  is  computed  from  the  eigenvalue  n.  If  this 
frequency  differs  sufficiently  from  that  of  the  previous  design,  the 
newly  calculated  frequency  is  used  to  recalculate  k,  and  the  analysis 
is  repeated.  (This  procedure  implies,  of  course,  that  [GA]  has  been 
determined  for  an  appropriate  range  of  k values,  as  discussed  above.) 


The  current  constraint  value,  g 
calculated. 


and  either 


13. 
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4.  NUMERICAL  EXAMPLES 

B'or  comparison  purposes,  three  separate  redesign  algorithms  were 
coded — two  based  on  equations  (6)  and  (7),  and  an  "energy-density" 
one  based  on  equations  (6) , (22) , and  (23) . The  first  two  differ 

in  the  selection  of  a and  in  the  treatment  of  side  constraints;  one 
(Segenreich)  is  based  on  the  algorithm  in  references  1 and  2,  while 
the  other  (Rizzi)  is  based  on  the  algorithm  in  references  6 and  7. 
These  are  described  in  subsection  2.3.  The  third  algorithm  ("energy- 
density")  is  described  in  subsection  2.5.  A fourth  algorithm  (ref. 
13),  based  on  the  method  of  feasible  directions,  was  also  included 
in  order  to  provide  some  comparisons  with  an  MP  method.  Program 
CONMIN,  described  in  reference  13,  can  use  analytically  computed 
gradients  or  can  calculate  gradients  by  finite  differencing.  Since 
analytical  gradients  were  already  being  calculated,  they  were  used 
for  CONMIN  as  well. 

4.1  Rectangular  Wing 

The  first  example  problem  involves  a simple  rectangular  wing 
structure  whose  dimensions  are  given  in  figure  1.  This  wing  was  first 
used  in  reference  14  and  has  since  been  treated  by  other  researchers. 

A very  simple  finite-element  model  was  created,  involving  two  cover 
sheets,  two  spar  webs,  one  rib,  and  four  spar  caps  in  each  of  three 
bays  in  the  structural  box.  The  spar  caps  were  represented  by  axial 
elements,  and  the  other  members  were  represented  by  in-plane  elements. 
There  are  12  design  variables,  whose  numbers  and  initial  values  are 
given  in  Table  1.  In  all  cases,  minimum-gage  constraints  of  one 
quarter  of  these  values  were  imposed.  No  weight  or  stiffness  is 
assigned  to  any  portion  of  the  wing  except  the  structural  box,  and 
the  initial  weight  is  88.45  kg.  This  is  the  initial  configuration 
for  all  of  the  constraints  considered.  All  of  the  computing  was 
performed  on  an  IBM  370/168  computer. 

For  the  displacement  constraint,  transverse  loads  of  445  N were 
applied  at  the  six  nodes  on  one  side  of  the  wing,  and  the  transverse 
displacement  at  the  tip  nodes  was  calculated  to  be  1.465  cm.  This 
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displacement  was  constrained  to  be  1.524  cm.  Iteration  histories  for 
the  various  algorithms  used  are  given  in  figure  2,  and  relative  CPU 
times  are  given  in  Table  2,  along  with  the  final  values  of  the  active 
design  variables,  the  weight,  and  the  constraint  parameter 
c = ^'^r  ^ des  ~ final  designs  are  essentially  iden- 

tical, with  the  only  active  variables  being  the  web  thicknesses  in 
bays  1 and  2.  The  Rizzi  and  "energy-density"  algorithms  took  less 
than  one  third  the  CPU  time  taken  by  CONMIN,  with  the  Segenreich 
algorithm  having  an  intermediate  time.  The  iteration  histories 
reveal  that  all  of  the  OC  algorithms  nearly  reach  the  optimum  weight 
within  ten  iterations.  However,  it  was  necessary  to  have  these 
algorithms  continue  in  order  to  ensure  that  all  of  the  passive 
variables  were  identified  and  the  constraint  satisfied.  For  all 
except  the  "energy-density”  algorithm,  the  intermediate  designs 
followed  the  constraint  boundary  rather  closely,  since  such  behavior 
is  enforced  at  each  step.  The  "energy-density"  algorithm  is  somewhat 
looser  in  this  respect  and  appears  to  proceed  in  two  stages  - a first 
stage  where  most  of  the  weight  is  removed,  and  a second  where  the 
constraint  is  satisfied.  Parameters  chosen  for  the  various  algorithms 
are  as  follows:  For  the  Segenreich  algorithm,  the  sequence  of  weight- 

reduction  factors  K.  was  (0.2,  0.1,  0.05,  0.025,  0.01,  0.005);  for 

^ ( 0 ) 

the  Rizzi  algorithm,  a = 0.90,  = 0.95;  for  the  "energy-density" 

algorithm,  e^^  = 0.5,  e^  = 2.0;  for  all  of  thes.  algorithms,  the  con- 
vergence parameter  e.  was  0.001;  and  for  CONMIN,  either  the  default 
values  or  the  values  recommended  in  reference  13  were  used  for  the 
various  parameters  required. 

For  a frequency  constraint,  the  fundamental  frequency  of  free 
vibration  was  constrained  to  be  68.0  rad/sec,  or  10.82  Hz.  This  is 
slightly  higher  than  the  calculated  frequency  of  67.16  rad/sec,  or 
10.69  Hz,  for  the  initial  design.  Since  both  the  inertial  and  stiff- 
ness properties  of  the  wing  are  linearly  proportional  to  the  design 
variables  and  there  is  no  nonactive  mass  or  stiffness,  constraining 
the  frequency  to  be  identical  to  that  of  the  initial  design  would 
result  in  a trivial  problem,  since  the  mass  and  stiffness  matrices 


[ 


-18- 


can  be  scaled  by  an  arbitrary  factor  without  affecting  the  frequency. 

The  optimum  woxild  therefore  bo  given  by  the  design  with  all  design 
variables  at  their  minimum  values.  Kith  a slightly  altered  frequency 
constraint,  the  optimum  design  should  have  at  least  one  active  design 
variable,  and  this  in  fact  is  the  case,  as  Cvin  bo  seen  in  Table  3. 

The  constraint  parameter  c is  here  ~ 1*  The  iteration 

histories  are  given  in  figure  3.  For  the  Segenreich  algorithm,  the 
table  of  weight-reduction  factors  was  (0.1,  0.05,  0.025,  0.01,  0.005, 
0.002);  for  the  Rizzi  algorit)un,  = 0.95,  = 1.0;  the  parameters 

for  the  "energy-density"  algorithms  and  for  CONMIN  were  identical  to 
those  used  for  the  displacement  constraint.  The  CPU-time  comparison 
shows  that  this  time  only  the  "energy-density"  algorithm  is  faster 
than  CONMIN,  although  the  differences  are  not  great.  The  "energy- 
density"  algorithm  actually  results  in  an  increase  in  weight  for  a 
few  iterations  before  it  converges,  and  the  Kizzi  algorithm  does  not 
display  the  same  rapid  approach  to  the  vicinity  of  the  optimum  weight 
as  it  did  for  the  displacement  constraint.  The  latter  phenomenon  is 
undoubtedly  a result  of  keeping  a constant , which  was  necessary  in 
order  to  achieve  convergence. 

For  the  flutter  constraint,  the  six  transverse  modes  of  the 
initial  design  were  used  to  define  generalized  coordinates,  and  sub- 
sonic generalized  aerodynamic  forces  were  calculated  with  the  doublet 
lattice  method  (the  program  described  in  reference  12  was  used) . A 
flutter  Mach  number  of  0.717  was  calculated  at  an  altitude  of  1,372  m, 
corresponding  to  a speed  of  240  m./sec.  This  flutter  p>oint  was  imp>osed 
as  the  constraint,  with  c = g.  CPU-time  comparisons  among  the  four 
algorithms  and  the  final  design  information  are  given  in  Table  4,  » nd 
the  iteration  histories  are  given  in  figure  4.  The  parameters  useu 
in  the  various  algorithms  are  the  same  as  those  used  for  the  displace- 
ment constraint,  with  two  exceptions  ~ c ~ 0.002  in  the  Rizzi  algorithm, 
and  0.1  in  the  "energy-density"  algorithm.  The  exponent  e^  had  to 

be  reduced  from  the  originally  selected  value  of  0.5  in  order  to  prevent 
divergence.  The  iteration  history  reflects  this  reduction,  in  that  the 
approach  to  the  optimum  weight  is  more  gradual  than  it  was  in  the 
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previous  examples.  As  a chock  on  the  accuracy  of  using  fixed  modes, 
the  fir\al  design  was  reanalyzed  for  flutter  with  new  natural  modes. 

This  new  analysis  gave  a flutter  speed  6%  lower  than  that  calculated 
with  the  original  modes.  Since  the  maximum  number  of  transverse 
modes  available,  six,  was  used  in  bot.h  cases,  it  was  not  possible  to 
see  whether  retaining  more  modes  initially  would  narrow  this  difference. 

<1.2  Rectangular  Wing  with  Increased  Degrees  of  Freedom  and  Internal 

Fuel 

In  an  effort  to  obtain  a slightly  more  complex  problem  with  more 
design  variables  active  at  the  optimum,  the  rectangular  wing  was 
remodeled  with  an  increased  number  of  elements,  as  indicated  in 
figure  5.  Each  bay  was  divided  into  two  with  a set  of  nodes  at  the 
midspan  of  the  bay.  New  ribs  were  not  added,  however,  so  only  the 
number  of  cover  sheets,  spar  webs,  and  spar  caps  was  doubled.  This 
resulted  in  a total  of  21  design  variables,  whose  numbering  and  initial 
values  are  given  in  Table  5.  In  addition,  nonstructural  mass  to  rep- 
resent internal  fuel  was  distributed  uniformly  within  the  structural 
box.  The  initial  weight  was  increased  in  this  manner  by  110.58  kg 
to  a total  of  199.13  kg.  Minimum-gage  constraints  of  25%  of  the 
initial  values  were  again  imposed  on  the  design  variables. 

The  initial  design  was  analyzed  for  flutter  with  generalized 
coordinates  defined  by  12  transverse  vibration  modes,  and  a flutter 
speed  of  231  m/sec  was  calculated  at  an  altitude  of  1,372  m,  correspond- 
ing to  a Mach  number  of  0.689.  As  before,  doublet-lattice  generalized 
aerodynamic  forces  were  calculated  with  the  program  described  in  refer- 
ence 12.  This  flutter  speed  was  imposed  as  the  behavioral  constraint, 
and  the  wing  was  optimized  with  the  Rizzi  algorithm,  with  = 0.90, 

= 0.95,  and  e = 0.005.  The  iteration  history  is  shown  in  figure  6, 
and  the  final  design  information  is  given  in  Table  6.  In  addition  to 
the  webs  that  were  active  before  at  the  optimum,  there  are  now  some 
cover  sheets  that  are  also  involved,  but  most  of  the  design  variables 
are  still  at  their  minimum  values.  The  final  design  was  also  reanalyzed 
for  flutter  with  new  modes,  and  the  recalculated  flutter  speed  is  9% 
lower  than  that  calculated  with  the  original  modes.  In  this  case,  also, 
all  of  the  modes  available  were  used  for  both  calculations. 
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Attompti!  to  obtairi  optimum  solutions  with  the  oilier  algorithms 
wore  not  succossrul,  and  at  tins  point  it  was  decided  to  use  another 
example  that  would  be  more  representative  of  a leal  design  and  that 
might  also  have  more  active  desiipi  vatiaLilcs.  This  is  described 
be  love. 

4. J  Delta  Wing 

The  next  example  structure  to  be  considered  is  a biconvex-airfoil 
delta  wing,  shown  in  planform  in  figure  7.  This  wing  is  example  2 of 
reference  15.  It  can  have  up  to  102  design  variables,  and  lower  totals 
can  be  created  through  linking,  so  problems  of  varying  complexity  in 
terms  of  the  optimization  task  can  be  treated  while  the  fvilly  modeled 
structure  is  used  for  analyses.  The  structural  planform  is  built  up 
from  triangular  in-plane  elements  separated  by  very  stiff  axial  members 
at  each  mode  that  model  the  core.  Mass  to  represent  internal  fuel  is 
distributed  in  the  core. 

The  structure  has  been  modeled  and  is  still  in  the  process  of 
being  analyzed,  so  there  arc  no  optimization  rv'sult.s  tJiat  can  be  given 
hero . 

5.  CONCLUDING  RKMARKS 

5.1  ConclusLonii  From  Results  to  Date 

From  the  results  wltlv  the  simple  model  of  the  rectangular  wing, 
it  can  be  seen  that  the  OC  algorithms  generally  perrormed  better  thatv 
the  M!'  algorithm  in  terms  of  relative  CFU  time,  with  the  Rizzi  algoriltim 
doing  the  best.  However,  the  need  for  the  user  to  «’xperimont  with 
these  parameters  in  any  new  appliciition  will  have  a great  influence 
on  tlve  utility  of  the  algorithm  being  used,  which  could  well  overshailow 
any  advant.njes  the  algorithm  might  have  in  computational  efficiency. 

For  the  final  comparisons,  therefore,  it  would  bo  desirable  to  liave 
recommi>ndpd  or  default  values  for  the  parameters  in  various  applications 
and  not  change  them  unless  it.  is  absolutely  necessary. 
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In  retrospect,  it  has  also  become  apparent  that  the  converycnice 
criteria  for  the  various  alyorithins  may  not  have  been  as  consistent 
as  they  should  have.  Consider,  for  example,  the  relationship  between 
the  convergence  criterion  for  the  Rizzi  algorithm  (cq.  (14))  and  t)iat 
for  the  "energy-density"  and  Segenreicli  algorithms  (eg.  (12)).  The 
Rizzi  and  Segenreich  criteria  can  be  related  through  the  recursion 
relation,  equation  (7).  This  can  be  manipulated  to  yield 
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If  is  the  convergence  parameter  for  the  Segenreich  algorithm  and 
that  for  the  Rizzi  algorithm,  then 


With  = 0.90  and  = 0.95,  after  25  iterations  a = 0.263,  so 

= 0.737e2*  the  other  hand,  if  = 1.0,  then  = 0.1e2- 

This  means  that,  for  equivalent  convergence  criteria,  in  these  two 
cases,  shovild  be  74%  or  10%  of  G2*  With  the  rectangular-wing 
example,  inconsistencies  in  these  criteria  affect  CPU  time  more  than 
the  final  answer,  since  the  optimum  appears  to  bo  located  in  a very 
shallow  depression  in  design  space.  It  is  hoped  that  the  delta-wing 
structure,  or  others  that  may  be  analyzed,  will  have  more  sharply 
defined  optima. 


5.2  Future  Work 

During  the  next  year,  comparisons  of  the  various  algorithms  will 
bo  continued  with  the  delta-\%'ing  example  described  above.  This 
example  will  also  be  used  to  test  the  accuracy  of  computiiig  gradients 
with  fixed  modes.  The  most  promising  algorithm  will  then  bo  extended 
to  treat  multiple  behavioral  constraints,  and  alternative  strategies 
for  handling  these  constraints  v/ill  be  evaluated. 
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TABLE  1.  DESIGN-VARIABLE  NUMBERING  AND  INITIAL  VALUES  FOR  THE  RECTANGULAR  WING 


TABLE  3.  RELATIVE  CPU  TIMES  AND  FINAL  DESIGN  INFORMATION  FOR  THE  OPTIMAL 
RECTANGULAR  WING  WITH  A FREQUENCY  CONSTRAINT. 
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TABLE  4.  RELATIVE  CPU  TIMES  AND  FINAL  DESIGN  INFORMATION  FOR  THE  OPTIMAL 
RECTANGULAR  WING  WITH  A FREQUENCY  CONSTRAINT. 
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DESIGN- VARIABLZ  NUMBERING  AND  INITIAL  VALUES  FOR  THE  RECTANGULAR  WING 


TABLL  6.  FINAL  DESIGN  INFORMATION  FOR  THE  RECTANGULAR  WING  WITH 
INCREASED  NUMBER  OF  ELEMENTS  AND  A FLUTTER  CONSTRAINT. 

Weight,  kg:  143.2 

Active  design  variables,  cm;  = 0.5469,  tg  = 0.3195,  tg  = 0.07341 

= 0.1569,  = 0.1301,  tj^^  = 0.04879 

-4 

Constraint:  -0.715x10 

Iterations:  24 
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Figure  7.-  Layout  of  delta-wing  structure, 
All  dimensions  are  in  cm. 
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Abstract 


It  is  observoJ  that  modern  optimal  design  of  structures  represents 
a confluence  of  two  streams  of  theoretical  development:  Matrlc  finite- 
element  approximation  on  the  digital  computer  — a technology  of  which 
Professor  Argyris  is  one  of  the  founders;  and  practical  application  of 
the  variational  calculus.  The  present  paper  addresses  optimization 
problems  wherein  complicated  constraints  involving  dynamic  aeroelastic 
behavior  are  prominent.  Search  procedures  based  on  optimality  criteria 
are  believed  to  offer  special  advantages  relative  to  such  problems. 

With  the  principal  constraint  formulated  in  terms  of  the  "V-g  method" 
of  flutter  analysis,  three  search  schemes  are  applied  to  the  minimum- 
weight  redesign  of  a particular  wing.  The  first  scheme  is  based  on  the 
method  of  feasible  directions  and  is  representative  of  mathematical- 
programming  methods.  The  other  two  are  deriv'ed  from  necessary  conditions 
for  a local  optimum  and  can  be  classed  as  optimality-criteria  schemes. 
Although  the  results  are  by  no  means  definitive,  they  do  suggest  that  a 
heuristic  redesign  algorithm  based  on  an  optimality  criterion  may  be  the 
best  candidate  for  Incorporation  in  a m.ore  general  design  procedure 
capable  of  treating  multiple  constraints  with  large  numbers  of  design 
variables. 

The  paper's  final  section  undertakes  to  shovj  how  optimality  criteria 
might  bo  constructed  when  the  aeroelastic  constraint  is  written  in  the 
time  domain.  Three  special  forms  of  the  aerodynamic  generalized  forces 
are  considered:  quasl-stcady , quasi-steady  with  dissipation  omitted, 
and  fully  unsteady.  The  resulting  criteria  for  the  first  and  last  cases 
are  based  on  an  unproven  hypothesis,  but  it  is  suggested  that  their 
simplicity  merits  a trial  application. 
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The  design  of  realistic,  efficient  structures  by  procedures  based 
on  the  mathematics  of  optimization  is  today  a routine  tool  in  the  manu- 
facture of  large  aircraft  and  other  devices  where  a combination  of  light 


weight  with  high  reliability  is  necessary  for  success.  Thif?  situation 
comes  about  in  consequence  of  the  merging  of  two  streams  of  analytical 
development.  The  first  of  these  involves  the  use  of  finite-element 
approximations  to  built-up  structural  arrangements,  together  with  matrix 

I 

i theory  and  the  digital  computer.  Although  hundreds  of  names  might  be 

i 

I mentioned,  no  one  has  contributed  more  to  the  foundations  or  to  the 

I 

I current  useful  state  of  finite-element  methods  than  Professor  Argyris. 

In  addition  to  numerous  research  papers,  his  several  books  (e.g., 

Argyris  [1])  and  his  valuable  series  of  articles  with  Professor  Kelsey 
In  Aircraft  Engineering  are  classics  of  the  field.  The  computer  program 
ASKA,  developed  with  his  colleagues  at  Institut  fur  statik  und  Dynamik 
der  Flugkonstruktionen  in  Stuttgart,  is  in  daily  use  throughout  the  world 
on  problems  which  include  the  class  discussed  in  this  paper. 

The  second  stream  derives  from  variational  calculus  and  the  concept 
of  extrema.  Based  on  the  construction  of  necessary  or  sufficient  conditions 
that  must  be  met  by  the  optimum  (usually,  minimum-weight)  design.  Its 
practical  aspect  consists  in  the  formulation  of  ever-morc-ef f Icicnt  search 
methods,  which  are  intended  to  bring  a trial  or  starting  configuration  to 
within  acceptable  convergence  of  the  goal  by  the  smallest  number  of 
steps  or  the  least  cycles  of  computer  operation.  An  excellent  summary, 
with  Innumerable  examples,  of  the  mathematics  of  optimization  is  contained 

in  Bryson  and  Ho  [2].  The  AIAA  Structures  Design  Lecture  by  Schmit  13) 

i 

reviews  the  rich  history  of  aeronautical  applications. 


The  present  research  addresses  a snail  corner  of  dynanic  structural 
opt inization:  one  where  constraints  relating  to  flutter  of  a wing  or  other 
djaiamic  aeroelastic  perfomance  r.ust  be  imposed  along  with  conditions  of  a 
more  conventional  nature,  such  as  those  relating  to  stress  under  load, 
deflection,  minimum  dimensions  of  structural  elements,  etc.  This  general 
topic  was  reviewed  recently  by  Stroud  [4]  — a survey  which  perhaps 
relieves  the  present  authors  of  responsibility  for  extensive  literature 
citations.  Special  recognition  should  be  given,  however,  to  the  paper  by 
Turner  [5] , wherein  optimization  with  rigorous  flutter  requirements  was 
first  formulated  in  a discrete,  finite-element  framework. 

The  focus  here  Is  on  a single  constraint  involving  aeroelastic 
stability.  Section  II  begins  with  a very  familiar  statement  of  the 
flutter  problem  for  a linear  system  with  a finite  number  of  degrees  of 
freedom  (cf.  Bisplinghoff  ct  al.  [6],  Sect.  9-5).  The  structure’s  motion 
Is  assumed  in  advance  to  be  a simple  harmonic  time  function,  and  through  the 
artificial  Introduction  of  energy  dissipation  one  seeks  the  actual  speed 
of  neutral  stability  for  flight  under  given  atmospheric  conditions. 

Because  flutter  calculation  is  so  time  consuming  in  cases  of  aero- 
nautical Interest,  there  is  here  a special  reason  for  identifying  search 
methods  that  require  this  step  as  Infrequently  as  possible.  It  is  the 
authors'  opinion  that  schemes  which  fall  under  the  heading  of  optimality 
criteria  offer  the  best  prospect.  Accordingly,  a relatively  simple  wing 
structure,  subjected  to  a single  constraint  on  its  flutter  performance. 

Is  analyzed  in  several  ways  and  the  computer  costs  are  compared.  The 
chosen  methods  range  from  woll-knov,-n  and  generally-available  routines, 
based  on  mathemat ical  programming,  to  a pure  criterion  approach  that 
is  believed  to  incorporate  some  new  features.  Wiile  the  latter  is 
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similar  in  form  to  those  employed  by  Siegel  [7]  and  Haftka  et  al . [8]  , 
it  alone  is  capable  of  converging  to  a true  local  optimum.  As  will  be 
seen,  this  permits  comparisons  of  efficiency  among  the  different  schemes 
that  involve  both  the  same  initial  design  and  essentially  the  same  final 
design. 

Section  III  contains  no  numerical  results  but  is  an  exploration  of  how 
aeroelastic  optimization  might  be  carried  out  in  circumstances  when  it  is 
undesirable  to  prescribe  simple  harmonic  motion.  Such  computations 
appear  more  feasible  today  because  of  investigations  line  those  of  Vepa 
[9]  and  Edwards  [10] , wherein  means  are  described  for  adapting  existing 
aerodj-namic  theory  to  the  unsteady  flows  produced  by  general  small  motions 
of  wings  or  bodies. 

The  aim  of  Sect.  Ill  is  to  produce  optimality  criteria  under  various 
approximations  to  the  aerod>Tiarcic  teinas  in  the  equations  of  notion. 

Although  no  such  schece  nay  be  regarded  as  proven  until  after  its  success- 
ful application  to  meaningful  designs,  nevertheless  these  proposals  are 
deemed  worthy  of  trial.  In  the  process  of  their  development,  the  concept 
of  the  adj.^-'  system  plays  a significant  role.  A curious  discovery'  is 
mentioned,  wnich  relates  to  this  adjoint  in  circumstances  where  "aero- 
dynamic memory"  must  be  accounted  for. 

II . C0^tPARIS0NS  OF  DIFFERENT  OPTIMIZATION  METHODS  ON  A VING  STRUCTURE 

WITH  FIXED  FLUTTER  SPEED 

In  his  lecture,  Schmit  [3]  has  discussed  the  different  philosophies 
underlying  optimization  schemes  derived  from  mathematical  programming  and 
from  optimality  criteria.  Basically,  the  mathematical-programming  algorithm 
make  use  of  information  from  the  current  design  and  calculate  a design 
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chanse  that  will  alter  the  objective  in  the  required  direction  without 
violating  any  of  the  constraints.  Optinal ity-criterion  algorithms  are 
derived  — often  purely  heurist ically  — from  the  necessary  conditions 
for  optimality,  or  an  approximation  thereto.  Either  approach  may  or 
nay  not  require  that  derivatives  of  the  constraints  be  calculated.  For 
a complicated  behavioral  constraint,  such  as  a constraint  on  flutter 
speed,  the  computation  of  derivatives  may  introduce  unnecessary  penalties 
In  terras  of  computation  tine.  The  purpose  of  this  section  is  to  introduce 
an  algorithm  that  does  not  require  derivatives  and  to  present  some  compari- 
sons with  other  algorithms  on  a relatively  simple  system  governed  by  a 
single  equality  constraint  on  the  flutter  speed.  Additionally,  the  new 
algorithm  will  be  seer  to  have  the  capability  of  converging  to  the  "exact" 
optimum,  so  that  the  comparisons  are  more  meaningful.  Although  virtually 
every  new  or  revised  scheme  that  has  appeared  has  been  compared  with  other 
schemes,  there  are  very  few  instances  (e.g.,  Haftka  et  al.  [8]  ) where  the 
comparisons  have  been  made  on  the  same  computer  system. 

2.1.  Statement  of  the  Flutter-Speed  Constraint 

Consider  now  a lifting  surface  whose  deformation  is  approximated  by 
superposition  of  a finite  number  of  free-vibration  modes,  v/hich  may  or 
may  not  be  normal  nodes.  (Typically,  they  will  be  the  normal  inodes  of 
an  initial  design,  which  for  the  sake  of  simplicity  will  be  retained  as 
primitive  modes  during  the  optimization.)  These  modes  define  generalized 
coordinates  i = 1,  ...,  n-  With  the  assumption  of  simple  harmonic 

motion  In  time  t , the  governing  equations  for  flutter  become 


M]  + (A]  - n [K] He}  = (0) 


)(C1 


(2.1) 
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Here  b'l  and  [k1  are  syrjnetric  nxn  generalized  mass  and  stiffness 
matrices,  respectively,  [a]  is  a matrix  of  oscillatory  generalized 
aerodynamic  forces,  and  C^(t)  - llie  matrix  [A]  is  a complex 

function  of  Mach  number  and  reduced  frequency.  In  the  "V-g"  method, 

(2.1)*is  solved  for  fixed  Mach  number  and  a number  of  reduced  frequencies 
to  give  a set  of  eigenvalues  Sl^,  and  eigenvectors  » i •=  1,  n. 

With  O + ,(j  ),  a corresponding  set  of  frequencies 

and  artificial  damping  parameters  can  be  calculated.  From  each 

frequency  and  the  reduced  frequency,  an  airspeed  V can  be  determined, 
and  the  roots  of  (2.1)  can  then  be  plotted  as  curves  of  V vs.  g for 
varying  values  of  reduced  frequency.  The  lowest  speed  at  which  a root 
makes  the  transition  from  negative  g,  denoting  stability,  to  positive  g, 
denoting  instability,  is  the  critical  flutter  speed.  A final  step  involves 
repeating  this  process  for  other  values  of  Mach  number  until  the  values  of  the 
critical  flutter  speed,  the  Mach  number,  and  the  sound  speed  at  the  chosen 
altitude  are  compatible. 

For  the  purposes  of  optimization,  there  are  several  ways  in  which 
the  constraint  may  be  imposed.  One  possibility  is  to  work  directly  with 
the  speed  itself;  this  has  been  done  successfully’  in  a number  of  instances 
( [ll],  [12]).  Another  possibility  that  offers  certain  advantages  is  to 
fix  the  speed  (and  therefore  the  Mach  number)  and  to  constrain  the  value 
of  g (cf.  Segenreich  and  McIntosh  [13]).  This  practice  will  be  followed 
here,  and  the  constraint  for  flutter  is  then  stated  simply  as 

g 1 0 (2.2) 


6 


for  the  critical  root  (or  node),  with  speed  and  Mach  nunber  fixed.  In 
the  exanples  that  follow,  this  constraint  will  be  the  only  behavioral 
constraint,  so  (2,2)  can  really  be  regarded  as  an  equality  constraint. 

2.2.  A Heuristic  Redesign  Algorithm  Based  on  an  Optinality  Criterion 

The  finite-element  representations  of  taany  structures  involve  elements 
whose  stiffness  and  Inertial  properties  depend  linearly  on  the  design 
variables,  which  are  conmonly  plate  thicknesses,  spar-cap  cross-sectional 
areas,  etc.  This  linear  dependence  will  therefore  be  assumed  here,  and 
the  weight  to  be  minimized  will  be  written  in  the  form 


J 

o 


N 


E 


"lb 


(2.3) 


There  may, in  fact, be  soma  mass  that  is  invariant  with  respect  to  the 
optimization,  but  it  is  not  necessary  to  Include  it  in  (2.3).  In  accor- 
dance with  Turner  [5]  and  Haftka  et  al. [8] , the  equations  of  motion 
(2.1),  written  for  the  critical  node,  are  prerr.ultiplied  by  a row  matrix 
of  Lagrange  multipliers  or  adjoint  variables,  and  the  real  part  of  this 
quantity  is  adjoined  to  to  give 


J({q},{?},  , t^)  = + Re  ( [q]  [B]{e})  (2. A) 

where 

[B]  = [M]  + [A]  - [K]  (2.5) 

Necessary  conditions  for  optimality  are  given  by  the  vanishing  of  variations 
of  J with  respect  to  n,  t^,  and  the  elements  of  (q)  and  . These 
yield,  respectively. 


NtfUU 


7 


(2.6) 

(2.7) 

(2.8) 

(2.9) 


Since  n » 1/oj  on  the  constraint  boundary  where  g ^ 0,  (2.6)  Is  equiv- 
alent to  the  vanishing  of  the  variation  of  J with  respect  to  the  flutter 
frequency  w , and  it  can  be  viewed  as  a relation  giving  the  flutter 
frequency  at  the  optimum  design.  The  original  constraint  equations  are 
reproduced  by  (2.8),  while  (2.9)  defines  the  adjoint  equations.  The 
optinulity  criterion  is  given  by  (2.7).  Under  the  aforementioned  assump- 
tion of  linear  dependence  of  the  inertial  and  stiffness  properties,  iM] 

and  [K]  can  be  written  an 
N 


[M] 

i“l 

ft 

IK) 

" t 
1=1 

and  (2.7)  becomes 


(2.10) 

(2.11) 


(2.12) 

The  left-hand  sides  of  (2.12)  resemble  energy  densities,  nnd  the  optimality 
criterion  is  therefore  seen  to  require  that  all  of  these  "energy  dctinltles" 
h.nvo  the  same  value.  If  mlnimnm-gage  constraints  on  the  design  v.arlables 
are  also  specified,  then  some  design  variables  may  become  passive — i.c., 
equal  to  their  minimum  allowable  va lues  — -“dur Ing  redesign.  If  thl.s  occurs, 
then  (2,1?)  must  Iiold  only  for  the  active  design  v.-ir  tables. 


a 


One  can  create  any  number  of  iterative  redesign  scheracs  based  on 
(2.12).  One  particularly  attractive  candidate  is  the  following: 


v+1 


'"i  ^i 


(e  ) le  e 

_JLL-  ' d + g'')  2 


(2.13) 

(2.14) 


av 


Here  v is  the  iteration  number,  and  e is  the  average  of  the  (e  ) 

av  Vi 

for  all  the  active  variables.  The  absolute  value  of  the  "energy-density" 
ratio  is  required,  since  the  expar.ent  e^  is  typically  less  than  unity 
and  the  ■’••X  either  positive  or  negative.  Note  that  as  the 

"energy-density"  ratios  approach  unity,  the  constraint  factor  (1  + g )**- 
serves  to  ensure  that  the  constraint  g = 0 will  be  satisfied.  The  form 
of  (2.14)  is  derived  from  two  assumptions: 

(a)  If  ld^)j|  > I t corresponding  design  variable 

should  be  increased,  and 

(b)  If  the  current  design  is  not  feasible  (g'^  >0),  all  design 
variables  should  be  increased. 

Convergence  of  this  iterative  formula  cannot  be  proven,  and  there  is  no 
guarantee  that  the  formula  will  he  capable  of  equalizing  not  only  the 
magnitudes  but  also  the  signs  of  the  d^)^  » which  is  a necessary  condition 
for  optimality  according  to  (2.12).  The  formula  (2.14)  is  similar  to  that 
used  by  Haftka  et  al.  [8]  , except  that  the  absolute  value  rather  than  the 
real  part  was  used  in  defining  the  (2.12).  In  effect,  this 

means  that  the  algorithm  used  in  [8]  attempts  to  satisfy  an  approximate 
rather  than  an  exact  optimality  criterion,  and  the  final  designs  obtained 
with  this  algorithm  did  not  coriespond  to  the  final  designs  obtained  v.'ith 
other  methods. 
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The  analysis  and  redesign  algorithm  proposed  herein  proceeds  in 
the  following  steps: 

(1)  The  current  value  of  the  constraint  is  found  by  solving  (2.8) 
and  obtaining  g,  w , and  (C)  for  the  critical  mode.  Since 
the  aerodynamic  matrix  [A]  will  vary  with  frequency,  it  is 
required  that  the  frequency  used  In  defining  [A]  coincide  with 
the  frequency  calculated  from  fl  . This  is  achieved  iteratively. 
The  frequency  from  the  previous  design  is  used  initially  to 
define  [A]  and  is  then  compared  with  the  frequency  calculated 
from  n . If  these  two  frequencies  are  not  in  agreement  to  within 
a specified  limit,  the  frequency  computed  from  is  used  to 
determine  a new  generalized  aerodynamic  matrix  [A]  and  the 
process  is  repeated. 

(2)  The  adjoint  equations  (2.9)  are  solved,  and  {q}  for  the 
critical  mode  is  obtained. 

(3)  The  densities  (e  ) '^  are  calculated  as  in  (2.12)  for  the 

i 

V 

current  active  set  of  design  variables,  and  e and  the 

" av 

V V 

C.  are  calculated.  The  C.  are  calculated  for  all  i . 
i X 

(4)  The  new  active  design  variables  ^ calculated.  If 

any  of  these  is  less  than  its  specified  minimum  value,  it  is 
set  to  that  miniraumvalue  and  is  relegated  to  the  passive  set. 

(5)  If  for  a passive  variable  > 1.0,  this  variable  is  reintro- 

duced to  the  active  set  and  steps  (3)  - (5)  arc  repeated  until 

the  active-passive  identities  are  stable.  Once  this  stability 

v+1 

is  achieved,  the  new  set  of  design  variables  t^  is  taken 

as  the  next  design. 


• 4 
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for  the  active 


(6)  Convergence  is  checked  by  testing  the 

set.  If  I 1 - C.'*  I < e for  all  of  the  tested,  the  new 

design  is  declared  final.  If  not,  the  algorithm  is  repeated. 

The  convergence  parameter  e is  set  by  the  user. 

Before  a numerical  example  is  given,  it  will  be  instructive  to  consider 
briefly  an  alternative  optimality  criterion  and  to  explore  its  relationship 
to  the  criterion  derived  above.  This  criterion  is  derived  from  variations 
of  the  function 


J(X.  t^)  - + Ag  (2.15) 

Vanishing  of  the  variations  of  J with  respect  to  the  design  variables 
yields  the  optimality  criterion 


3t, 


1 

\ 


(2.16) 


Once  again,  this  equality  holds  for  those  design  variables  that  are  active 
at  the  optimum.  Following  Segenrelch  and  McIntosh  [13]  , one  can  write  for 
the  derivatives 


. [(r/  -«2R/  -gl‘)(2sR3-f2l3+  ^ I^) 

- (2R^  - 2g  I3  + ^ R^)  + g R2^)]  /D  (2.17) 

where 

3 3 

D *=  (2g  R^  + 2I3  + I^)  I3  + (2R3  - 2g  I3  + ^ R^)  R3,  (2.18) 

Rj^  = Re  ( iqj  [M^]  iO  ),  (2.19) 

Rj^  = Re  ( i^i  [K^]  U)  ),  (2.20) 
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(2.21) 


= Re  ( iqj  [K]{r}  ), 

= R ( iqj  [l^jU)  ),  (2.22) 

i i 

Ij  j Ij  , and  are  the  corresponding  imaginary  parts,  b is  a ref- 

erence length,  and  k = wb/V  is  the  reduced  frequency.  At  first  glance, 
it  would  seem  that  (2.16)  bears  little  resemblance  to  its  counterpart  (2.12). 
However,  by  making  use  of  (2.6),  which  is  not  obtained  from  the  formulation 
immediately  above,  it  can  be  shovm  that  at  the  optimum 


2R3  - 2g  I3  + 


bcj' 


(2.23) 


Making  use  of  this  relation  in  (2.18)  and  (2.17)  produces  a vastly  simplified 
expression  for  , which  when  inserted  into  (2.16)  gives 


1 , 2_  i ^ i . ^ i . 

— (a.  Rj  - R2  + g I2  ) = 
i 


(2.2A) 


_ 2 

Since  Q = (1  + jg)/a)  , (2.12)  can  be  rewritten  with  the  definitions  (2.19)- 

(2.22)  to  give 

12  11  1 2 

— (m  Rj  - R^  + g I2  ) = - w (2.25) 

Since  the  quantities  , etc.  involve  triple  matrix  products  with  {q} 

and  (C)  , and  since  the  normalization  of  (q)  and  is  arbitrary, 

(2.24)  and  (2.25)  are  in  fact  equivalent  expressions.  It  must  be  emphasized, 

3 o 

however,  that  the  simplified  expression  for  v—  is  valid  only  at  the 
optimum,  so  that  an  optimization  scheme  that  requires  derivative  calculations 
must  use  the  full  expression  (2.17), 


12 


r 


2.3.  A Numerical  Example 

In  order  to  illustrate  sonie  of  the  ideas  discussed  above,  consider  the 
rectangular  wing  whose  diraensions  are  given  in  Figure  1.  This  wing  was 
created  by  Rudisill  and  Bhatia  [11]  and  has  since  been  considered  by 
Segenreich  at)d  McIntosh  [13],  among  others.  The  structural  box  has  three 
bays.  In  each  bay  there  are  two  cover  sheets,  two  spar  webs,  and  one  rib, 
all  modeled  by  in-plane  rectangular  elements,  and  four  spar  caps,  modeled 

by  axial  elements.  The  design  variables  are  numbered  as  shown  in  Table  1. 

2 

The  initial  design  values  were:  tj^  - = 12.90  cm  , t^  - t^-  » 0.2032  cm., 

t^  - tj2  0.1016  cm.  Tne  weight  of  the  initial  design  was  88.45  kg. 
Minimum-gage  constraints  were  imposed  at  one  quarter  of  the  Initial  values. 
For  the  flutter  calculations,  the  first  six  transverse  f ree-vlbration  modes 
of  the  initial  design  were  used  to  define  generalized  coordinates.  The 
doublet-lattice  method  [14]  was  used  to  calculate  generalized  aerodynamic 
forces,  and  a flutter  Mach  number  of  0.717  was  calculated  at  an  altitude 
of  1,372  meters.  This  flutter  point  was  imposed  as  the  primary  behavioral 
constraint.  Optimal  designs  were  calculated  with  three  methods  — the 
"energy-density"  method  described  in  3.2,  the  method  used  by  Segenreich  and 
McIntosh  [13],  and  a method  developed  by  Vanderplaats  [15].  The  method  in 
[13]  is  derived  from  an  algorithm  due  to  Kiusalaas  [16],  which  is  based  on 
an  optimality  criterion  of  the  form  (2.16)  and  requires  calculation  of  the 
derivatives  at  each  Iteration.  The  method  in  [15]  is  based  on  a 

mathematical-programming  procedure  known  as  the  method  of  feasible  directions 
[17]. 

The  final  designs  from  all  three  methods  were  essentially  identical. 

The  only  active  design  variable  was  t^,  the  thickness  of  the  front  and 
rear  spar  webs  in  bay  1.  The  final  designs,  constraint  values,  and  relative 
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computer  CPU  times  are  given  in  Table  2.  All  computations  were  performed 
on  an  IBM  370/168  computer  at  Stanford  University.  Both  Siegel  [7]  and 
Haftka  et  al.  [8]  used  e^  = 0.5,  e^  = 2.0  in  their  computations.  With 

these  values,  the  algorithm  of  Sect,  2.2  diverged,  and  it  was  necessary  to 
reduce  e^^  to  0.15  in  (2.14)  before  convergence  could  be  demonstrated. 

Before  the  flutter  constraint  was  considered,  the  same  structure  was 
optimized  with  a displacement  and  a natural-frequency  constraint.  For  these 
problems,  it  was  possible  to  leave  e^  at  0.5.  All  three  methods  produced 
almost  identical  final  designs  and  CPU-time  comparisons  similar  to  those 
in  Table  1.  ^ - 

The  final  design  was  reanalyzed  for  flutter  with  new  (normal)  modes 
in  order  to  ensure  that  the  root  that  was  constrained  was  still  the  critical 
root.  This  was  found  to  be  the  case;  the  critical  branches  of  the  Initial 
and  final  designs  are  compared  in  Fig.  2.  Using  the  normal  modes  of  the 
initial  design  as  primitive  modes  during  redesign  is  seen  to  cause  an 
error  of  5.5%  in  the  estimation  of  the  flutter  speed  in  comparison  with 
the  flutter  speed  calculated  with  normal  modes  of  the  final  design. 

In  more  complex  practical  applications  of  the  method,  especially  when  the 
initial  design  is  far  from  optimal,  one  should  occasionally  re-calculate 
the  normal  modes.  They  arc  then  to  be  used  as  the  basis  for  subsequent 
Iterations,  and  slight  mismatches  of  the  final  flutter  speed  can  thus  be 
avoided. 
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III.  OPTIMUITY  CRITERIA  BASED  ON  AEROEL/vSTIC  EQUATIONS  IK  THE  TIME 
ANP  LAPLACE-TRANSFORiM  DOMAINS 

Since  optimality  criteria  are  closely  related  to  variational  prin- 
ciples from  which  the  associated  equations  of  notion  can  be  derived,  the 
first  stage  of  this  investigation  looks  at  these  principles  and  at  the 
adjoint  systems  which  are  an  inherent  feature  when  nonconservative  effects 
are  present.  The  starting  point,  as  in  Sect.  II,  is  a reference  design 
whose  behavior  can  be  adequately  represented  by  the  generalized  coordinates 
C^Ct)  , 1 = 1,  ...  n,  of  some  subset  of  its  natural  inodes  of  free  vibration. 

The  generalized  masses  and  stiffnesses  are  described  by  symmetric  nxn 
matrices  [M]  and  [K] , respectively.  At  the  beginning  of  optimization,  these 
matrices  may  be  diagonal.  But,  as  the  search  proceeds,  the  generalized 
coordinates  and  associated  eigenfunctions  are  left  unchanged.  Therefore 

and  [KCrij^)]  become  full  matrices,  functions  of  the  N added  or  sub- 
tracted masses  that  constitute  the  vector  of  design  variables.  The 

dependence  on  mj^  inay,  in  one  or  both  instances,  be  linear  (cf.  Turner  [5]), 
but  this  is  not  necessary  to  the  development. 

For  the  analysis  of  homogeneous  aeroelastic  phenomena  like  flutter 
stability,  the  equations  of  motion  governing  the  column  matrix 
contain  aerodynamic  generalized  forces,  linearly  dependent  on  this  matrix. 

Two  cases  are  considered  in  what  follows. 

3. 1.  Quasi-Steady  Aerodynamics 

In  this  approximation  (cf.  Sect.  5-6  of  Bisplinghoff  et  al.  [6] 
and  Pines  [18],  among  many  other  sources  of  information)  there  is  negligible 
"mem.ory,"  so  that  the  airloads  are  algebraically  related  to  ^ and 


^The  "apparent  mass"  effect  might  also  be  included  through  a 1^1  term,  but 
this  refinement  v.’ould  contribute  nothing  of  significance. 
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tM]  {O  + lK]{0 


([D]{C 


{|}  + IE]{U 


) ^ ° 


(3.1) 


Here  [D]  and  [E]  are  nonsymnetric  catrices  of  real  constants  for  a given 
flight  condition.  It  proves  desirable  to  extract  from  these  terms  a 
quantity  with  dimensions  of  force,  containing  the  aeroelastic  eigenvalue 
This  eigenvalue  may  be  the  actual  flight  dynamic  pressure  or  perhaps 
Mach  number,  etc. [cf . the  parameter  used  by  Welsshaar  ([19] , Chapter  6) 

in  panel  flutter  optimization ).  Equations  (3. 1)  then  become 


[M]  {C)  + [K]{5>  - S cA;  Id)  CO  + Ie){5} 


where  [d]  and  [e]  would  normally  be  dimensionless.  S is  a reference 
area  such  as  that  of  the  wing  planform  projection;  reference  wing  chordlengtu 
c and  flight  speed  V are  employed  to  cancel  the  dimensions  of  the  time 
derivative. 

Let  q^  (t) , j = > n,  be  generalized  coordinates  for  the  adjoint 

system.  The  equations  governiiig  the  column  matrix  of  the  q^  are 


tM]{q}  + [K]  {q}  - Q S ( - c/V  Id]"^  {q}  + [e]^  {q} 


).0 


(3.3) 


Here  superscript  T denotes  the  transpose  of  a square  matrix.  One  way 
of  deriving  (3.3)  is  to  require  stationarlty,  with  respect  to  variations 
in  the  elements  of  (C)  and  (q),  of  a generalized  Hamilton’s  principle 
which  might  be  written 

H = ^ - Iqj  + IqJ  [KK^}  + Q^S  ]lql  [d]{C} 


- Q S iqj  leK?}  >dt 


(3.4) 
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The  open  brackets  mean  a row  m.atrix.  In  the  usual  way,  "initial"  con- 
ditions are  prescribed  at  the  time  limits  t^  and  to  prevent  nonzero 
terms  from  arising  during  integration  by  parts.  It  is  remarked  that  H — 
like  the  optimality  criteria  based  on  it,  to  be  proposed  below  — is  not 
unique.  Thus  the  third  term  in  the  integrand  could  have  its  sign  changed 
and  the  time  differentiation  transferred  from  iqj  to  . 

3.2.  Linear  Aerodynamics  with  "Memory." 

When  analyzing  stability,  experience  shows  that  one  must  usually 
account  for  the  unsteady  influences  of  the  past  history  of  a wing  or  body’s 
motion.  They  are  due  to  the  presence  of  a vortex  wake  shed  as  the  result 
of  prior  changes  in  the  circulation  "bound"  to  a wing,  to  the  finite  speed 
of  sound  propagation  in  the  gas,  or  to  both  of  these  effects.  One  must 
then  replace  the  aerodynamic  terms  in  (3. 1)  by  a convolution 


[A]  * iO  = 


(A(t  - t)]  (C  (t)}  dx 


(3.5) 


In  this  form,  the  motion  is  assumed  to  have  begun  at  t = 0.  The  elements 
AfjCt)  of  [a]  are  indiclal  time  functions,  giving  the  generalized  force 
exerted  on  one  degree  of  freedom  i due  to  impulsive  motion  in  another  j. 
In  most  cases  these  elements  are  known,  without  approximation,  only  in 
terms  of  their  Laplace  transforms  A^^^  (s)  . 

Replacing  [a]  by  a dimensionless  equivalent  [a]  as  before,  one  is 


led  to  the  fully  unsteady  generalization  of  (3.2); 


[Mice)  + K{C}  - Q^S 


la(t-x)l  {^(t)}  dx 


(3.6) 


Clearly,  (3.6)  can  be  derived  from  stationarity,  with  respect  to  the 


elements  of  iqj , of 
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< - IqJ  [mKC)  + IqJ  [kKC) 


/' 


~ iq(t)j  / [a(t-T)l  {C(t)}  dt  > dt 


(3.7) 


It  is  of  interest  to  examine  the  adjoint  equation,  the  one  governing 
{q(t)}.  To  this  end  and  before  a ceaningful  variation  of  H with  respect 
to  can  be  carried  out,  the  order  of  integrations  in  (3.7)  must  be 
Inverted  and  the  variables  t and  x interchanged.  Since  the  result  should 
not  depend  on  the  choice  of  t^  and  t^*  it  seems  convenient  to  replace 
these  quantities  with  the  starting  time  tj  = 0 and  with  ^ • After 

these  choices  have  been  made,  simple  manipulations  lead  to 

[M]  {q}  + [K]  {q}  - Q S / (a(T  - t)] {q(T) }dT  = 0 (3.8) 


I 


Albeit  there  are  other  ways  of  writing  the  adjoint  equation,  (3.8) 
has  a certain  intrigue.  If  one  attempts  to  give  physical  meaning  to  the 
aerodynamic  term,  it  would  seem  to  imply  that  the  state  at  instant  t is 
affected  by  events  subsequent  in  time.  Such  speculation  is  apparently'  not 
fruitful.  As  in  other  nonself-adjoint  systems,  there  is  often  a great 
deal  of  artificiality  in  the  adjoint  problem,  and  this  feature  is  here 
exacerbated  by  the  presence  of  a convolution  in  (3.6).  The  convolution 
theorem  of  Laplace  transformation,  incidentally,  does  not  apply  to  the 
integral  in  (3.8). 

In  an  effort  to  gain  insight,  an  elementary  problem  with  a single 
degree  of  freedom  C(t)  was  constructed  from  (3.7)  and  (3.8).  The  jingle 
indicial  function  was  taken  as  a simple  lag. 


A(t)  = a e 
o 


- Ot 


(3.9) 


Even  in  this  case,  the  solution  of  (3.6)  and  initial  conditions 
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I 

i 

I 


+ M:j  C - a 


f 


- o (t  - t)  . , V , 

e ' 5 (t)  dT 


C(0)  = ; C(0)  - 0, 


(3. 10a) 
(3.10b) 


Is  not  particularly  convenient  to  vrite  out.  The  special  choice 
2 

a = Ka  = Mil)  o causes  a pole-zero  cancellation,  however,  yielding 
o 


the  elementary  result 


C(t)  (1)  ot/2  . , ^ , N 

— — = — e sin  (it)  t +ij,'), 

C 0)  n 

o n 


(3.11a) 


with  (i)  = ii) 

n 


/l  - (o/2(d)‘ 


tp  = tan  ^ (2ti)  /a) 
n 


(3.11b) 

(3.11c) 


The  corresponding  adjoint  q(t),  also  with 


q(0)  = q^;  q(0)  * 0 


(3.12) 


! 

K 

t 


cannot  be  calculated  by  Laplace  transformation.  Merely  assuming  it  in 
the  form  of  (3.11a),  however,  yields 

q(t)/q^  = 5(t)/5^  (3.13) 

This  is  physically  reasonable  and  not  unexpected,  because  the  system  and 
its  adjoint  should  have  Che  same  eigenvalues. 

3.3.  The  Search  for  Optimality  Criteria. 

As  an  aid  to  deriving  criteria  that  might  enable  efficient  search 
routines  while  aeroelasClc  constraints  are  enforced,  the  ideas  in  Sect.  5b 
of  Plant's  review  [20]  furnished  an  excellent  lead.  This  presentation  is, 
accordingly,  broken  into  three  parts,  the  first  dealing  just  with  a discrete 
version  of  the  same  kind  of  nonconservative  system  treated  by  that  author. 


i 

i 
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3.3a.  Minimam-UeiKht ^ tructure  w h Prescribed  Flutter  S tability 
Boundary  of  the  Froguency-Merslng  Type . 

In  the  quasi-steady  Eq.  (3.2),  specialize  further  by  assuming 


[dl  = 0 


(3. 1/i) 


This  elitiination  of  dissipative  mechanisms  ensures  that  the  eigenvalues  of 


+ [kH^}  - Q^S  [e]{C>  = 0 


(3.15) 


remain  purely  imaginary  as  parameter  is  increased  from  zero.  Because 

of  the  asymmetry  of  [e] , however,  there  is  some  at  which  a pair 

becomes  equal  and  produces  "instability  by  merging  of  frequencies"  (cf. 
Pines  [18]  ). 

In  a form  slightly  different  from  that  in  Sect.  II,  a typical  optimi- 


zation problem  may  be  stated  as  follows: 


Minimize  ^ ^ , 

k=l 

with  {5}  a solution  of  (3.15)  such  that 


Q«c  ~ % 


(3.16) 


(3.17) 


and  (for  example) 


1 ^ = 1.2.  .. . N . 

o 


(3.18) 


Observe  that  the  can  be  inferred  from  specified  minimum  dimensions 

*^o 

of  structural  members;  if  represents  the  reduction  of  total  weight  from 
an  initial  design,  they  will  be  negative  numbers.  Equation  (3.15),  its 
adjoint,  and  the  generalized  Hamiltonian  Integrand  (cf.  3.^0  arc  rewritten  w V 
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the  substitutions 

s *=  - 0, 

^i(t)  = q e®'' 

and  q^Ct) 

St 

e ; 

(3.19) 

1 - Q [M]  + [K]  - 

a)  = 

0 

(3.20) 

|-f>  [m]  + [K]  - 

{q}  = 

0 

(3.21) 

h = iqj  [K]  (T) 

+ fi  iqj 

[M]  (U 

- Q S w [c]{i'} 

9 

(3.22) 

[Note  that,  although  H is  nonlinear  In  q and  5,  it  is  expected  that 
the  eigenvectors  will  be  real  up  to  Therefore,  there  is  no  question 

of  taking  real  parts  in  (3.22).] 

Plaut  [20]  points  out  that  the  characteristic  equation  of  (3.20)  or 
(3.21)  has  the  form 

Q„)  = 0 (3.23) 

He  reasons  that,  with  two  frequencies  merged  at  f>  = and  Q = 
it  should  be  true  that 


d Q 
^00 

“dsT 


0 


He  suggests  promultlplylng  (3.20)  by  Iqj  and  solving  for  : 


-niq]  [M]{cj  4-  tqj  [K]{c} 
S iqj  [e] 


(3.24) 


(3.25) 


Applying  (3.24)  to  (3.25)  then  produces  a relation 
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- 0 , 


(3.26) 


iql  [M]  (f.) 


which  mvist  hold  at  the  critical  condition. 

The  manner  of  const  runt  1 n>»  an  optimality  criterion  Is  to  adjoin  h 
to  the  performance  index  with  a hagrange  multiplier  X , 


) “ £ 

' k"l 

+ xf[q]  tK(nj^)){?}  + iUq]  lM(mj^)]{f,>  - Q^S  IqHeKu]  (3.27) 


The  variations  (or  partial  derivatives)  of  J with  respect 
of  (t)  and  {q}  produce  (3.21)  and  (3.20),  respectively, 
the  variation  with  respect  to  0 is  assured  by  (3.26). 
on  the  design  variables  m^^  yield 


to  the  elements 
The  vanishing  of 
Finally,  variations 


(3.28) 


Note  that  the  equal  it  for  all  k " 1,  ...»  of  the  quantities  in 
parentheses  in  (3.28)  furnishes  the  desired  criterion.  Just  as  in  Sect.  II 
the.  "energy  ratios"  are.  employed  to  adjust  the  individual  m.isses  during 
progress  toward  an  optimum,  so  these  quantities  should  be  able  to  do  the  job 
here  without  the  need  for  differentiations  with  respect  to  eigenvalues. 

That  such  criteria  do  load  to  meaningful  designs  is  implied  by  the  simple 
examples  given  In  Plant  120]., 

3.3b.  More  General  Quasi-Steady  Airloads 

The  physical  and  adjoint  equations  of  motion  in  this  case  are  the 
full  (3.2)  and  (3.3).  To  ho  specific,  let  it  be  required  that,  at  a 
given  0^,  the  least  stable  of  the  neroelastlc  eigenvalues  has  a given 
degree  of  stability  — a condition  that  might  bo  expressed 
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Re  {s.}  < - s , for  all  1 = 1,2 
1—0 


(3.29) 


Here  is  a positive  real  constant.  Under  assumption  (3.19),  the  real 

or  complex-conjugate  pairs  of  s^  are  eigenvalues  of 


< s [M]  + [K]  - Q S 


[d]  + [e]  >{(}  - 0 


(3.30) 


and  also  of  the  corresponding  equations  from  (3.3).  Subject  to  these 
constraints  and  to  a prescribed  Q^,  the  optimization  problem  then 
consists  of  (3.16),  (3.18)  and  (3.29). 

By  analogy  with  (3.27),  it  is  now  hypothesized  that  an  appropriate 
formulation  requires  the  vanishing  of  variations  with  respect  to  q^. 


and  s of 


J ({?},  Cq},  s,  \ ^ + Re  <x(iqj[K]{C}  - lqj(M]{C} 

' ^ k= 1 ' 

- j iqj  [d]  {?}  + Q^S  iqj[e]{C}>  (3. 


(3.31) 


Here,  of  course,  X , the  eigenvectors  and  s^^  arc  all  complex  numbers. 
The  line  of  reasoning  whereby  the  real  part  is  taken  in  (3.31)  follows 
directly  that  given  first  by  Turner  [5]. 

As  before,  two  additional  relations  (beside  the  physical  and  adjoint 
equations)  are  obtained  from  variations  of  (3.31). 


- 1 - Ro<x(cqj[|;^]{0  - > 


(3.32) 


for  all  k ^ 1 N.  Also 


Re  <x|2s  lqj[H](f,}+  Q„S  ~ iqj  [d]{C>J  > ° 0 


(3.33) 
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The  solution  procedure  that  is  proposed  to  be  tried  starts  from 
calculating,  at  a given  design  step,  all  eigenvalues  and  vectors  of  (3.30) 
and  its  adjoint.  In  general  (3.32)  and  (3.33)  will  not  be  satisfied.  For 
a particular  eigenvalue  s (the  least  stable),  they  amount  to  N + 1 real 
equations  involving  N + 2 variables:  the  N X™  (a). 

Using  the  known  for  the  given  design,  one  might  select  a single  one 

of  (3.32)  plus  (3.33)  and  calculate  A . With  A available,  the  remaining 
members  of  (3.32)  can  be  used  to  adjust  the  various  mj^  up  or  down,  aiming 
for  a design  in  which  the  quantities  within  the  parentheses  of  (3.32)  are 
equal  for  all  k. 

One  difficulty  remains:  that  of  how  to  drive  the.  real  part  of  the 
least  stable  of  the  s^  tow'ard  s^  at  each  step.  This  might  be  done  by 
numerically  differentiating  the  stability  determinant  of  (3.30)  with 
respect  to  the  "most  active"  of  the  design  variables  tDj^.  This  mj^  could 
be  altered  by  just  the  amount  needed  to  bring  Re  to  Then  the 

aforementioned  optimality  criterion  from  (3.32)  provides  a w’ay  of  correcting 
all  other  rij^.  (Note  that  "passive"  m^,  which  already  coincide  with  their 
minima  prescribed  by  (3.18),  are  omitted.) 

Obviously,  the  changes  in  the  other  m^^  (besides  the  "most  active") 
are  going  to  throw  off  the  system  stability.  Nevertheless,  these  changes 
become  progressively  smaller  as  the  optimization  progresses,  so  choosing 
the  particular  index  k whose  mass  seems  to  be  varying  most  rapidly  may 
furnish  a sufficiently  powerful  control  on  Re  {s^}. 
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3.3c.  Fully  Unsteady  Airloads 

Surprisingly,  this  case  appears  to  be  a fairly  straightf orvard 
extension  of  the  treatment  in  Sect.  3,3b.  The  key  step  seems  to  be  to 
assume  that  the  motion  has  gone  on  for  a long  time  prior  to  t = 0. 

In  this  case,  (3.6)  and  (3.8)  must  be  modified  to 


[M]{C>  + IK]{C}  - Q^S 


I la(t-T 

•/.CO  OD 

L 


)]  {C(t)}  dr 


(3.34) 


[M]  (q)  + [K]{q}  - Q S / [a(T-t)]  {q(T)}  dx  = 0 


(3.35) 


Leaving  aside  certain  questions  of  existence,  one  again  chooses  the 
homogeneous  solutions  (3.19).  This  substitution  into  (3.34),  with  the 


variable  change  x^  = t - x,  then  leads  to  the  algebraic  system 


[K]  + s^  [M]  - Q S [a  (s) 


0. 


(3.36) 


where  the  a^^ (s)  are  the  Laplace  transforms  of  the  aerodynamic  indicial 
functions  — exactly  the  quantities  which  the  Edwards  I 10]  investigation 
shows  how  to  calculate. 

Inserting  (3.19)  into  (3.34)  gives  rise  to  the  following  manipulations 
of  the  aerodynamic  integral: 


00 

/- 

■if. 


la(x-t)]"  {q( 


sx. 


) (q) 


^s(x-t)  dxj{q} 


ta(x2)]'  dx2  I (q)  e®"  = [a  (-s)]"^  {q}  e®*" 


(3.37) 


In  (3.37)  X2  = X - .t,  and  the  elements  of  the  final  square  matrix  are 
simply  the  Laplace  transforms  of  » '^ith  the  sign  of  s reversed. 
Substituting  (3.37)  into  (3.35),  one  obtains 
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(3.38) 


|[K]  + [M]  - QJS  [a  = 0 

The  steps  in  the  proposed-optimality— criterion  approach  parallel 
rather  closely  those  starting  from  (3.31).  The  expressions  containing 
matrices  [d] , [e]  and  their  tranposes  must,  of  course,  be  adjusted  by 
the  insertion  of  ta(s)]  and  [a  (-s)]  , as  appropriate.  The  quantity 

whose  variations  yield  the  needed  data  is 

J = + Re  < A |iqj  [K]  {?}-  iqj[M]{T}  - Iq J U(s)  ] (oj  > 

k=l  (3.39) 

The  principal  complication  over  the  quasi-steady  case  seems  to  be  that 

(3.33)  is  replaced  by  a more  elaborate  formula  containing  derivatives  of 

the  elements  of  fa(s)].  Although  these  derivatives  might  have  to  be 

obtained  numerically  in  the  "exact"  formulation,  one  suspects  that  their 

computation  will  be  greatly  simplified  by  adopting  the  Pad6  approximants 

of  Vepa  [9] . 
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IV.  CONCLUDING  REMARKS 


Although  tho  exanple  considered  in  Sect.  V is  a very  simple  one,  the 
results  are  encouraging  enough  to  support  additional  investigations  with 
the  new  algorithm  described  in  Sect.  2.2.  It  will,  of  course,  be  necessary 
to  apply  all  of  the  methods  employed  in  Sect.  2.3  to  more  realistic  problems 
before  any  quantitative  assessments  can  be  given  concerning  relative 
efficiency.  In  particular,  it  should  be  noted  that  the  weight  reduction 
achieved  — some  74%  of  the  initial  weight  — is  not  truly  representative 
of  what  would  occur  in  actual  practice.  Ultimately  the  methods  need  to  be 
tested  for  multiple  behavioral  constraints.  The  mathematical-programming 
procedure  CONMIN  [15]  is  already  capable  of  treating  such  problems,  and 
th«^  votV.  of  Segenreich  et  al.  [21]  and  Rizzi  [22]  will  provide  interesting 
possibilities  for  introducing  a multiple-constraint  capability  into  the 
algorithm  of  Sect.  2.2. 

Raving  thus  demonstrated,  by  means  of  a familiar  example,  the  attrac- 
tiveness of  optimality  criteria  in  connection  with  an  analytically-complicated 
design  constraint  like  flutter,  the  paper  proceeds  into  a more  speculative 
area.  It  cannot  be  too  strongly  emphasized  that  expressions  like  (3.31) 
and  (3.39)  are  hypothesized.  The  similarity  of  criteria  derived  from  them 
to  proven  counterparts  like  (3.27),  however,  lends  credence  to  the  proposal 
that  these  schemes  be  accorded" a fair  trial.  One  may  cite  the  interesting 
and  successful  work  of  Siegel  [7]  as  evidence  that  complete  mathematical 
rigor  is  not  always  necessary  in  a procedure  for  design  optimization. 
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FIGURE  CAPTIONS 


Figure  1.  Layout  of  rectangular  wing.  All  dimensions  are  in  era. 


Figure  2.  Behavior  of  critical  modes  in  V-g  plane  for  initial  and  final 
designs.  M = 0.717,  altitude  1,372  m,  6 modes. 


